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1.  Introduction 

The  present  work  focuses  on  electromechanical  behavior  of  di¬ 
electric  crystalline  solids.  A  dielectric  is  an  insulating  material  that 
exhibits  polarization  in  the  presence  of  an  electric  field.  Electrome¬ 
chanical  behaviors  of  interest  include  piezoelectric,  pyroelectric,  and 
ferroelectric  effects.  Piezoelectricity,  in  a  general  sense,  refers  to  the 
coupling  between  electric  field  or  polarization  and  stress  or  defor¬ 
mation.  In  continuum  theories,  piezoelectricity  of  first  order  is  at¬ 
tributed  to  the  particular  choice  of  free  energy  functional  for  the 
body  that  may  depend,  for  example,  on  the  product  of  the  strain 
and  the  polarization.  Such  first-order  piezoelectric  effects  can  only 
occur  in  crystal  classes  that  lack  a  center  of  symmetry  [43].  Second- 
order  electromechanical  effects  can  arise  in  non-conductors  of  all 
crystal  classes  as  a  result  of  quadratic  influences  of  the  electric  field. 
This  phenomenon  is  often  called  electrostriction  [19,20].  Pyroelectric 
crystals  exhibit  surface  charges  when  uniformly  heated  or  cooled; 
such  crystals  feature  energetic  coupling  between  polarization  and 
temperature.  The  pyroelectric  effect  is  revealed  by  heating  a  crystal 
to  induce  a  change  in  its  polarization.  Ferroelectric  crystals  comprise 
a  subset  of  pyroelectric  crystals,  the  former  exhibiting  spontaneous 
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polarity  that  can  be  reversed  by  an  applied  electric  field.  Ferroelec¬ 
tric  crystals  may  exhibit  a  transition  temperature,  called  the  Curie 
point,  above  which  they  are  not  spontaneously  polarized;  the  loss 
of  polarity  may  accompany  a  polymorphic  phase  transition  from  a 
non-centrosymmetric  to  a  centrosymmetric  structure. 

Many  theories  of  geometrically  and  materially  non-linear  elec¬ 
tromechanics  of  dielectric  solids  have  appeared  in  historic  and 
more  recent  literature.  Stratton  [57]  and  Landau  and  Lifshitz  [43] 
presented  formulations  encompassing  both  electrostatics  and  elec¬ 
trodynamics,  though  large  deformations  of  the  material  were  not 
thoroughly  considered.  Devonshire  [20]  developed  a  continuum 
thermodynamic  theory  of  ferroelectric  crystals  accounting  for  mate¬ 
rial  non-linearity,  e.g.  a  higher  than  quadratic  dependence  of  the  free 
energy  upon  polarization,  but  not  geometric  non-linearity.  Toupin 
[59],  Eringen  [27],  and  Tiersten  [58]  formulated  theories  of  elastic 
dielectric  bodies  subjected  to  arbitrarily  large  deformations.  Tiersten 
[58]  also  considered  thermal  effects  and  material  inertia.  Mindlin 
[52]  developed  frameworks  accounting  for  spatial  gradients  of  me¬ 
chanical  strain  and  polarization  and  demonstrated  correlation  be¬ 
tween  higher-order  continuum  theory  and  discrete  lattice  dynamics 
in  the  limit  of  long  wavelength  behavior.  Chowdhury  et  al.  [10]  for¬ 
mulated  a  non-linear  theory  for  thermoelastic  dielectrics  with  effects 
of  polarization  gradients.  Geometrically  non-linear  theories  of  elec¬ 
tromechanics  were  also  posited  by  Maugin  [48,49],  Hadjigeorgiou 
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et  al.  [30],  Dorfmann  and  Ogden  [22],  McMeeking  et  al.  [50],  and 
Vu  and  Steinmann  [60].  In  some  dielectrics,  reasonably  large  strains 
are  feasible  via  domain  switching  [64],  necessitating  the  use  of 
geometrically  non-linear  theory.  Large  deformations  are  also  at¬ 
tained  when  pressures  are  significant  enough  to  suppress  fracture, 
for  example  in  shock  physics  experiments  [49]  wherein  confined 
piezoelectric  ceramic  crystals  such  as  silicon  carbide  [51]  are 
encountered. 

Lattice  defects  are  known  to  strongly  affect  electromechanical 
behavior  of  dielectric  solids  and  hence  performance  of  engineering 
devices  fabricated  from  such  materials.  For  example,  dislocations 
accommodate  misfit  strains  between  dielectric  thin  films  and  sub¬ 
strates  in  electronic  devices  [56].  In  ionic  crystals  such  as  corun¬ 
dum  and  sodium  chloride,  consideration  of  charge  distributions  and 
effects  of  electric  fields  becomes  necessary  for  describing  disloca¬ 
tion  motion  and  dislocation  reactions  [34,40,41].  Polarized  domains 
and  domain  walls  affect  the  hysteresis  behavior  and  performance 
of  ferroelectric-based  actuator  systems  [64].  Dislocations  [5]  and 
vacancies  [8]  have  been  observed  in  quartz,  a  piezoelectric  crystal 
used  frequently  for  pressure  transducers  and  resonators.  Regard¬ 
ing  the  latter  type  of  defect,  vacancies  are  observed  in  a  number 
of  dielectric  materials  and  are  a  focus  of  the  theory  developed  in 
the  present  work.  Mobile  vacancies,  in  conjunction  with  climbing 
dislocations,  can  dominate  creep  deformation,  often  preferentially 
to  glide-controlled  inelasticity  at  high  temperatures  [9,61].  In  ionic 
crystals,  vacancies  typically  carry  an  electric  charge  [18,37].  Such 
charged  vacancies  notably  influence  dielectric  properties  and  elec¬ 
trical  loss  characteristics  of  capacitors,  oscillators,  and  tunable  fil¬ 
ters  [19],  for  example  those  comprised  of  perovskite  ceramic  crystals 
such  as  barium  titanate  and  strontium  titanate  [16,53].  Electroni¬ 
cally  active  Si  and  C  vacancies  are  important  in  silicon  carbide  [3],  a 
wide  band-gap  semiconductor. 

Some  clarification  of  terminology  used  for  defects  is  now  in  order. 
From  an  atomistic  perspective,  a  vacancy  is  regarded  as  an  empty 
atomic  site  in  the  crystal  structure,  i.e.  a  missing  atom.  From  the 
perspective  of  discrete  defects  in  elastic  continua  [28],  a  vacancy  can 
be  treated  as  a  sphere  on  the  order  of  the  atomic  volume  inserted  into 
a  slightly  larger  hole  in  the  solid.  The  boundary  of  the  hole  is  then 
pulled  into  rigid  contact  with  that  of  the  sphere,  resulting  in  a  radial 
displacement  field  in  the  solid  that  decays  with  distance  from  the 
defect.  A  void,  on  the  other  hand,  consists  of  multiple  missing  atoms 
and  is  typically  treated  in  continuum  frameworks  as  a  spherical  hole, 
free  of  traction  along  its  surface,  with  dimensions  significantly  larger 
than  that  of  a  vacancy. 

In  the  present  work,  a  theory  for  dielectric  solids  undergoing  po¬ 
tentially  large  volume  changes  from  mobile  vacancies  is  developed, 
with  defects  capable  of  carrying  an  electric  charge.  Large  deforma¬ 
tions  resulting  from  generation  and  motion  voids  or  vacancies  are 
important  when  considering  large  defect  concentrations  in  the  vicin¬ 
ity  of  grain  boundary  depletion  layers  [55],  and  in  general  when 
considering  coalescence  leading  to  fracture  [35].  A  finite  deforma¬ 
tion  theory  of  void  nucleation  and  growth,  also  incorporating  devia- 
toric  dislocation  plasticity,  was  developed  by  Bammann  and  Aifantis 
[2]  for  ductile  metallic  crystals,  but  diffusion  and  electromechanics 
were  not  considered.  Finite  volumetric  deformations  arising  from 
point  defects,  including  vacancies  and  interstitials,  were  addressed 
in  differential-geometric  treatments  of  Kroner  [42]  and  Clayton  et  al. 
[13],  but  electromechanical  effects  were  not  considered.  Many  and 
Rakavy  [46]  considered  charge  transport  resulting  from  diffusion  of 
various  carriers,  which  may  include  point  defects,  but  did  not  con¬ 
sider  mechanical  deformation.  Electrochemical  potentials  for  diffu¬ 
sion  of  charged  atomic  species  (i.e.  self-diffusion  or  mass  transport), 
interstitials,  and  vacancies  have  been  the  focus  of  a  number  of  stud¬ 
ies  [37,40].  The  electromechanical  behavior  of  dielectric  solids  con¬ 
taining  mobile  vacancies  was  more  recently  addressed  by  Xiao  and 


co-workers  [62,63]  and  Clayton  et  al.  [14,15].  Related  material  mod¬ 
els  have  accounted  for  surface  diffusion  [14,15,32]  and  dislocations 
[14].  Unlike  the  current  work,  none  accounted  for  finite  volumetric 
deformations  associated  with  vacancy  content. 

In  the  present  work,  electromechanics  of  dielectrics  in  the  context 
of  quasi-static  electric  fields  is  considered.  In  the  quasi-electrostatic 
approximation  [49,58],  finite  material  velocities  are  considered,  but 
electrodynamic  terms  in  Maxwell’s  equations  are  not.  Material  ve¬ 
locities  are  restricted  to  remain  small  compared  to  the  speed  of  light. 
Free  charge  densities  apart  from  charges  associated  explicitly  with 
vacancies,  for  example  free  electrons  and  electron  holes,  are  assumed 
quasi-static  within  the  dielectric  body.  This  is  in  contrast  to  an  elec¬ 
trical  conductor  with  free  electrons  whose  motion  may  lead  to  sub¬ 
stantial  current  flow,  requiring  a  formal  electrodynamic  description. 
The  vacancy  flux,  however,  can  be  interpreted  as  an  effective  current 
whose  divergence  is  proportional  to  a  rate  of  contribution  to  the  to¬ 
tal  free  charge  density  from  vacancies,  in  which  case  the  dielectric 
with  charged  defects  can  be  formally  classified  as  a  semiconductor. 
Neither  magnetic  effects  nor  mass  transport  (i.e.,  self-diffusion  of 
bulk  atoms  or  interstitials)  is  considered  here. 

The  remainder  of  this  paper  is  organized  as  follows.  First,  gov¬ 
erning  relationships  for  geometrically  non-linear  electrostatics— 
Maxwell’s  equations,  momentum  and  energy  balances,  and  the 
dissipation  inequality  in  the  context  of  finite  deformations— are  re¬ 
viewed.  Such  a  review  is  necessary  to  accompany  derivations  that 
follow.  A  geometrically  non-linear  theory  for  elastic  dielectrics,  first 
in  the  absence  of  defects,  is  presented.  Constitutive  relations  for 
elastic  dielectric  solids  emerge,  following  consideration  of  the  dis¬ 
sipation  inequality  [10,17,24,30,50].  Then,  following  similar  mathe¬ 
matical  and  physical  arguments,  a  geometrically  non-linear  theory 
of  dielectric  crystals  containing  mobile  vacancies  is  developed.  The 
vacancies  may  support  an  electric  charge  as  often  occurs  in  ionic 
solids  or  may  be  electrically  neutral  as  a  special  case  more  appli¬ 
cable,  for  example,  in  monatomic  crystals.  The  theory  presented 
here  extends  previous  work  [14,15]  to  large  deformations  and  to 
situations  wherein  Maxwell’s  stress  is  non-negligible,  though  sur¬ 
face  diffusion  and  surface  growth  considered  previously  are  not 
considered  here.  Rather,  particular  attention  is  directed  towards 
the  diffusion  equation  for  the  bulk  vacancy  flux.  Lastly,  effects  of 
large  elastic  deformations  and  large  vacancy  concentrations  on 
diffusion  of  vacancies  are  examined  for  a  slab  subjected  to  biaxial 
lattice  strain,  for  example  a  film  with  residual  stresses  due  to  lattice 
mismatch  [14,53]. 

The  following  notation  is  used.  Vector  and  tensor  quantities  are 
written  in  bold  type,  while  scalars  and  individual  components  are 
written  in  italics.  The  index  notation  follows  the  Einstein  summation 
convention,  distinguishing  between  covariant  (subscript)  and  con- 
travariant  (superscript)  components.  Current  configuration  indices 
are  written  in  lower  case  Latin,  reference  configuration  indices  in  up¬ 
per  case  Latin,  and  intermediate  configuration  indices  are  denoted  by 
Greek  symbols.  Juxtaposition  implies  summation  over  two  repeated 
adjacent  indices  (e.g.  (AB)ab  =AacBcb ).  The  dot  product  of  vectors  is 
represented  by  the  symbol  •  (e.g.  a»b  =  aagabbb ,  with  gab  compo¬ 
nents  of  a  metric).  Angled  brackets  denote  a  dual  (scalar)  product 
(e.g.  (A,  B)  =  tr( AB)  =AabBba  and  (a,b)  =  otaba).  The  colon  denotes  con¬ 
traction  over  repeated  pairs  of  indices  (e.g.  A  :  B  =  tr( ArB)  =  AabBab 
and  (C  :  A)ab  =  CabcdAcd ).  The  symbol  <g>  represents  the  tensor  prod¬ 
uct  (e.g.,  (a®  b)ab  =  aabb).  The  symbol  o  represents  the  composition, 
such  that  for  two  functions  /  and  g,  (f  og)(X)  =/(g(X)).  Superposed 
-1,7,  and  •  denote  inverse,  transpose,  and  material  time  deriva¬ 
tive  operations,  respectively.  Subscripted  commas  denote  partial  co¬ 
ordinate  differentiation.  Indices  in  parentheses  are  symmetric,  that 
is  2 A^  =  Aab  +  Aba,  while  indices  in  brackets  are  anti-symmetric, 
i.e.,  2 A^  =  Aab  -  Aba.  An  extensive  list  of  symbols  is  provided  in 
Appendix  A. 
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2.  Geometrically  non-linear  electromechanics 

Field  variables  associated  with  finite  deformations  and  electrome¬ 
chanical  effects  are  introduced  in  Section  2.1,  followed  by  governing 
equations  of  electrostatics  in  Section  2.2.  In  Sections  2.3,  2.4  and  2.5, 
equations  of  momentum  and  energy  conservation  and  the  dissipa¬ 
tion  inequality  for  dielectric  solids  subjected  to  finite  deformations 
and  electromechanical  loadings  are  presented. 

2.1.  Kinematics  and  electromechanical  quantities 


/s 


e 


Let  x=cp(X,  t)  denote  the  smooth  and  invertible  motion  of  a  body, 
with  x  denoting  spatial  coordinates  and  X  denoting  reference  coordi¬ 
nates.  The  deformed  body  B  is  parameterized  by  spatial  coordinates, 
while  the  body  Bo  in  the  reference  configuration  is  parameterized 
by  reference  coordinates.  The  local  deformation  gradient  is  the  map 
from  the  reference  to  current  tangent  space,  TxBq  TXB  [47] 


_  _  dxa  _A  _a  dxa 

F  =  Tcpx  =  — —  ga  <g>  G,  Fjt  =  — — — , 

x  dXA  A  dXA 


(1) 


where  ga  and  are  bases  referred  to  current  and  reference  coordi¬ 
nates,  respectively.  A  differential  line  element  dX  e  TXB 0  is  deformed 
to  dx  =  F  dX  e  TXB.  Mass  conservation  provides  the  requirements 


/  Podv  =  Pdv  ->•  p0  =  pj, 
Jv  Jv 


(2) 


where  reference  and  current  mass  densities  are  p0  and  p,  reference 
and  current  volume  elements  are  dV  and  dv ,  and  the  Jacobian  ]  = 
dv/dV  =  det(F“  X/g/G  >  0,  with  g  =  det (gab)  and  G  =  det(G^).  Spatial 
and  referential  metric  tensors  satisfy  gab  =  ga  #g b  and  gab  =  Ga  *GB. 
Christoffel  symbols  are 


2  r  bc  —  g  (Scd,b  +  Sbd,c  Sbc.d  )> 

2  T  =  Gad  (  GCd,b  +  Gbd,c  ~  Gbc,d  )•  ( 3 ) 

Components  of  the  covariant  derivative  of  a  spatial  vector  field  a  e  TB 

g  g 

are  then  vbaa=aab+r  bcac.  Analogous  formulae  apply  for  the  covari- 
g  ’ 

ant  derivative  Vb  taken  with  respect  to  reference  coordinates,  and 
covariant  derivatives  of  higher- rank  tensors  follow  conventional  def¬ 
initions  [13].  Specifically  let  v(x,t)  denote  the  spatial  velocity  field. 
The  spatial  velocity  gradient  is 


Fig.  1.  Continuum  quantities  for  deforming  dielectric  body  in  presence  of  electric 
field. 


Polarization  density,  or  simply  the  polarization,  is  the  vector  p  e 
TXB  (dimensions  of  charge  per  unit  area)  that  is  often  collinear  with 
the  relative  shift  of  electron  clouds  and/or  ions  comprising  a  di¬ 
electric,  or  with  the  orientation  of  a  permanent  electric  dipole  [49]. 
Polarization  vanishes  in  vacuum  and  in  conductors.  Electric  displace¬ 
ment  d  e  TB  is  related  to  electric  field  and  polarization  by  [57] 

d  =  e0e  +  p,  (6) 

where  the  dimensional  constant  e0(8.854(10)-12  C2  N-1  rrr2)  is  the 
permittivity  of  free  space.  Since  polarization  vanishes  in  vacuum,  (6) 
reduces  to  d  =  e0e  in  free  space.  The  surface  free  charge  density  a  is 
defined  on  surface  s  of  a  body  with  unit  normal  n  as 


ff=([d],n>.  (7) 

In  (7),  [d]  =  d+  -  d_  is  the  jump  in  d  at  the  interface  of  s  and  the 
medium  external  to  the  body,  which  could  be  vacuum  or  another 
body.  Here,  +  and  -  denote  respective  limiting  values  of  a  quantity 
at  locations  outside  and  inside  the  body  as  s  is  approached  from 
the  corresponding  side,  and  n  is  directed  from  the  -  side  to  the  + 
side  (directed  from  inside  to  outside).  A  deforming  dielectric  body 
in  three-dimensional  Euclidean  space  E3  with  corresponding  elec¬ 
tromechanical  quantities  is  shown  in  Fig.  1. 

2.2.  Maxwell’s  equations  of  electrostatics 

Maxwell’s  equations  of  electrostatics  consist  of  the  following  in¬ 
tegral  relations  in  the  spatial  frame: 


L%  =  V*  =  vafb  +  f  lacvc  =  FaAF~2A.  (4) 

Furthermore,  let  EAB = ( \  )(Qb  -  Gab  )  denote  components  of  the  sym¬ 
metric  Lagrangian  strain  tensor,  where  CAb  =  FaAgabFbB- It:  follows  that 
)  =JDaa  and  EAB  =  F^DabFbB,  where  Dab  =  L(ab). 

The  spatial  electric  field  e  e  TB  describes  the  Lorentz  force  f 
associated  with  a  point  charge  of  magnitude  q:  f =qe.  The  electric  field 
may  permeate  the  dielectric  body,  surrounding  media,  and  vacuum. 
The  electric  field  vanishes  inside  conductors  in  the  absence  of  current 
flow. 

The  free  charge  density  per  unit  spatial  volume  is  defined  by 

p  =  ^fiW(,(0  =  Efi(,WO, 

i  i 

with  the  number  of  charge  carriers  per  unit  volume  of  charge 
qh)=ez(z),  where  e  is  the  charge  of  an  electron  (1.602(1 0)_1 9  C)  andz(z) 
is  the  integer  valence  of  each  member  of  charge  carrier  population 
i.  For  excess  electrons,  z=  -1,  while  for  holes  or  missing  electrons, 
z=+l.  The  charge  density  p  vanishes  in  pure  vacuum  (i.e.,  a  vacuum 
containing  no  free  electrons)  and  within  neutral  conductors. 


I  e*dx  =  0,  f(d  ,n)  ds  =  f  pdv.  (8) 

Jc  Js  Jv 

The  first  of  Maxwell’s  equations  states  that  the  line  integral  of  the 
electric  field  along  an  arbitrary  closed  curve  c  vanishes.  From  Stokes’s 
theorem,  the  first  of  (8)  can  be  written 

^  ea  dxa  =  jT  sabcvaebnc  ds  =  0.  (9) 

A  vector  field  in  a  simply  connected  domain  whose  skew  covariant 
derivative  (i.e.  curl)  vanishes  can  be  represented  as  the  gradient  of  a 
scalar  potential.  Here,  the  scalar  potential  </)  is  called  the  electrostatic 
potential  and  is  continuous  throughout  space  [27].  The  local  form  of 
(9)  is1 

£abcVaeb  =  0  o  e[ba]  =  -tfciba]  =  0  o  h  =  ~&,b-  (10) 


1  Though  curvilinear  coordinates  are  used  throughout  the  paper  for  generality, 
Cartesian  indices  are  implied  for  all  integral  equations  involving  vector  and  tensor 
fields.  Alternatively,  following  [59],  all  quantities  entering  a  vector-valued  integrand 
can  be  parallel  transported  to  a  single  point  using  the  shifter,  and  the  resultant 
integral  evaluated  at  that  point. 
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In  the  domain  of  (9),  (^-continuity  of  e  is  assumed;  finite  jumps  in 
the  normal  component  of  e  are  admitted  across  surfaces  of  discon¬ 
tinuity,  but  £abclebJnc  =  0  across  such  surfaces  [27].  The  second  of 
Maxwell’s  equations  establishes  a  conservation  law  between  surface 
and  volumetric  charges.  Applying  the  divergence  theorem  to  the  left 
side  of  the  second  of  (8)  and  localizing  the  result 

Vada  =  p.  (11) 

In  arriving  at  (11),  (^-continuity  of  d  is  assumed  within  the  domain 
of  integration;  the  normal  component  of  a  jump  in  d  across  a  surface 
is  given  by  (7).  Multiplying  (11)  by  <J>,  integrating  over  volume  v 
enclosed  by  surface  s,  and  applying  the  divergence  theorem  with  (7) 
and  (10)  yields 

/  d *edv=  /  p(j)dv-\-  /  co^ds,  (12) 

Jv  Jv  Js 

where  co  =  -(d-,n)  measures  the  contribution  from  the  inside  of  s 
to  the  free  charge  density  in  (7). 

Spatial  relations  (8)-(12)  can  be  mapped  to  their  referential  coun¬ 
terparts  as  follows  [22,60].  In  particular,  the  first  of  (8)  becomes 


the  introduction  of  an  electromechanical  body  force.  In  electrostat¬ 
ics,  this  force,  per  unit  current  volume,  is 

b  =  pve  +  pe,  ba  =  pbvbea  +  pea,  (20) 

where  the  first  term  on  the  right  side  is  attributed  to  short  range  in¬ 
teraction  of  the  polarization  with  the  electric  field  [58,59],  and  the 
second  term  arises  from  Lorentz  forces  attributed  to  non-vanishing 
free  charges.  From  (6),  ( 1 0)  and  (11),  and  assuming  sufficient  smooth¬ 
ness  of  the  electric  field  and  polarization,  the  force  b  can  be  ex¬ 
pressed  as  the  divergence  of  a  rank  two  contravariant  tensor  known 
as  the  Maxwell  stress  x 

T  =  e®d -  y(e-e)g_1 

=  e®p  +  £0e®e- y(e-e)g_1,  vbzab  =  ba.  (21) 

With  the  inclusion  of  electrostatic  body  force  b,  the  global  balance 
of  linear  momentum  is 

f  pva  dv=  f  tads  +  f  (ba  +  ba)dv  =  f  (vbaab  +  ba  +  ba)dv,  (22) 

Jv  Js  Jv  Jv 


je>dx  =  j  e*FdX  =  J  Fre*dX  =  J  E-dX  =  0,  (13) 

where  C  is  a  closed  reference  curve  and  E=FTe  e  TB 0,  with  EA=F^ea. 
Relations  analogous  to  (10)  then  emerge 


£^bcWaEb  =  0  E[ba]  =  -&[,ba]  =  0  =  -<PB,  (14) 

where  the  potential  <£(X,  t)  =  <J>  o  cp.  From  Nanson’s  formula 
jf<D-,N)  dS  =  fvPodV,  (15) 


with  ta=aabnb  the  mechanical  traction,  a  the  Cauchy  stress,  and  b  the 
mechanical  body  force  per  unit  current  volume.  In  (22),  the  Cauchy 
stress  is  assumed  differentiable  within  v.  The  local  form  of  (22)  is 

Vbzab  +  ba  =  pva ,  (23) 

where  x  =  a  +  x  is  the  total  stress  tensor,  i.e.,  the  sum  of  the  Cauchy 
and  Maxwell  stresses.  Because  the  electric  field  and  polarization  may 
exhibit  jump  discontinuities  across  s,  the  boundary  conditions  are 

Ta  =  r°  -  [Ta>„,  (24) 


with  D  =  JF-1d  e  TBo,  p0  =  pj ,  and  5  the  surface  enclosing  reference 
volume  V.  The  analog  of  (11)  is  then  obtained  directly  from  (15)  and 
the  divergence  theorem 

VaD*  =  f'o-  (16) 

Finally,  multiplying  both  sides  of  (16)  by  &,  integrating  over  V ,  and 
using  (14)  gives 

[  D*EdV=  [  pQ&dV  -  [  <P(D_,N)  dS.  (17) 

Jv  Jv  Js 

No  natural  definition  exists  for  the  reference  analog  of  spatial  polar¬ 
ization  p.  One  obvious  assumption  is 

P  =  FTp,  (18) 


where  Ta  are  components  of  the  net  applied  traction  [59],  which  can 
be  assigned  as  Ta  =  t+a  [58]. 

The  body  force  b  likewise  enters  the  balance  of  moment  of  mo¬ 
mentum,  along  with  an  additional  moment  attributed  to  interaction 
between  the  polarization  and  electric  field  [58,59] 

Jv  zabc*bpvc  dv=  Jv  eabcxbbc  dv  +  j  sa6cxY'  ds 

+  /  eabcxbbcdv+  /  eabcecpbdv.  (25) 

Jv  Jv 

Application  of  the  divergence  theorem  and  Reynolds  transport  the¬ 
orem  along  with  (21)  leads  to 

f  zQbcxb(pvc  -bc  -vdTcd)dv  =  (  £abc{<jcb  +  pvcvb  +  tcb)dv.  (26) 
Jv  Jv 


leaving  the  following  referential  version  of  (6): 

D  =yC-1(s0E  +  P),  (19) 

where  the  symmetric  deformation  tensor  C  is  defined  immediately 
following  (4). 

Some  authors  partition  the  electric  field  or  electrostatic  potential 
into  contributions  from  various  internal  and  external  sources  [27,59]; 
in  the  present  work,  a  single  electric  field  and  potential  are  used  in 
the  governing  equations,  following  [58]. 

2.3.  Momentum  conservation 

Electromechanical  interactions  induce  modifications  to  the  bal¬ 
ance  laws  of  linear  and  angular  momentum  of  classical  continuum 
mechanics.  Such  interactions  can  be  addressed  straightforwardly  via 


The  left  side  of  (26)  vanishes  by  linear  momentum  balance  (23), 
and  the  second  term  in  the  integrand  on  the  right  vanishes  by  the 
symmetry  of  v  ®  v.  The  local  balance  of  angular  momentum  is,  from 
the  remaining  terms  in  (26), 

zab  =  z(ab\  erlai’l  =  zlbal  =  el  bpa|.  (27) 

meaning  that  the  total  stress  x  is  symmetric,  but  the  Cauchy  stress  <r 
need  not  be.  In  non-polarized  media,  =o  such  that  the  classical 
balance  of  angular  momentum  =  0  applies.  However,  if  such 
materials  support  an  electric  field  and  a  non-vanishing  free  charge 
density,  the  balance  of  linear  momentum  (23)  will  be  affected  by 
(20).  In  pure  vacuum,  the  Cauchy  stress  vanishes,  and  the  balance 

g  ,  yv 

of  linear  momentum  reduces  to  Vbt  =  ba  =  0,  satisfied  identically 
since  the  polarization  and  free  charge  density  vanish  by  definition 
in  pure  vacuum. 
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Various  definitions  have  been  postulated  for  the  mechanical 
stress  and  Maxwell  stress,  as  discussed  by  Eringen  [27].  The  defini¬ 
tion  used  in  the  present  work  for  the  Cauchy  stress  corresponds  to 
the  local  stress  of  Toupin  [59]  and  is  the  transpose  of  the  mechanical 
stress  of  Tiersten  [58].  The  definition  used  in  the  present  work  for 
the  Maxwell  stress  (21)  is  consistent  with  that  of  Toupin  [59]  and 
is  the  transpose  of  the  Maxwell  stress  of  Tiersten  [58].  Notice  that 
in  the  present  work,  the  index  of  the  traction  vector  corresponds  to 
the  first  index  of  the  corresponding  stress  tensor  (e.g.  ta  =  aabnb ), 
following  the  notation  scheme  of  Toupin  [59],  but  opposite  to 
notations  of  Eringen  [27],  Tiersten  [58],  and  McMeeking  et  al.  [50]. 

2.4.  Energy  conservation 

Various  methods  have  been  set  forth  to  account  for  energy  con¬ 
servation  in  dielectrics  [10,14,15,20,27,48-50,58-60].  Variational 
principles  can  provide  insight  into  field  equations  and  boundary 
conditions  for  static  and  non-dissipative  processes  [27,59],  while 
rate  forms  of  the  energy  balance  are  useful  for  situations  involving 
dynamics  and  dissipation  [14,15,50,58].  Prescribed  here  is  a  global 
balance  between  rates  of  internal  energy,  kinetic  energy,  external 
work,  and  thermochemical  heating 

^-(<f+jn=#+ x  (28) 

at 

with  the  contributions  from  kinetic  energy  Jf  and  extrinsic  thermo¬ 
chemical  energy  J  given  by 


where  application  of  the  chain  rule  produces  the  identity  ^  a  = 

4>jbVavb  -  e  •  Various  forms  have  been  suggested  for  Q  [27,49,58], 
and  some  authors  subtract  part  or  all  of  the  contribution  of  Q  from 
the  rate  of  total  internal  energy  rather  than  incorporate  it  in  the  rate 
of  external  power.  Identical  quantities  can  be  represented  in  a  vast 
number  of  ways  via  manipulation  of  Maxwell’s  equations  and  use 
of  vector  identities  and  theorems  of  Gauss  and  Stokes.  Recall  also 
that  in  the  present  Section,  along  the  lines  of  previous  theories  for 
dielectric  media  in  the  quasi-electrostatic  approximation  [10,50,58], 
purely  mechanical  dynamic  effects  are  considered  (i.e.  finite  veloc¬ 
ity  v),  but  electrodynamics  are  not  (i.e.,  fluxes  of  free  electrons/holes 
are  absent). 

Substituting  (29)-(32)  into  (28)  and  converting  all  surface  inte¬ 
grals  to  volume  integrals  using  (12)  and  the  divergence  theorem,  the 
global  balance  of  energy  becomes 

[  pedv+  f  pvavadv-  f  eapadi/|  +  f  ead  dv+ ^  f  eaeavbvbdv 
_Jv  Jv  Jv  J  Jv  2  Jv 

=  | J  prdv-  j  vbqb  dv+  J{vb(Jabva+(TabVbVa)dv 

+  J (ba+ba)vadv J+  J  ead  dy+y  J  (j) ^ v bvb  dv 

+  j{ea+^ya-4)fbVaVb+^,a^hyh+^aVbVb)dadv.  (33) 

The  integrand  of  the  last  term  on  the  right  side  of  (33)  vanishes 
identically.  Terms  in  braces  vanish  identically  in  vacuum,  and  terms 
not  in  braces  cancel.  Collecting  terms  in  braces  and  localizing  gives 


jf=  f  (p/2)vvdv,  J  =  (  prdv  -  S{ q  ,n)ds,  (29) 

Jv  Jv  Js 

where  a  dielectric  body  occupying  spatial  volume  v  with  oriented 
surface  element  nds  is  considered.  Field  variables  are  assumed  to 
possess  sufficient  smoothness  within  v  to  enable  use  of  local  forms 
of  Maxwell’s  equations  and  momentum  balances  as  well  as  the  di¬ 
vergence  theorem.  However,  electric  field  and  polarization  may  ex¬ 
hibit  finite  jumps  across  s.  In  the  present  thermodynamic  analysis, 
the  body  is  treated  as  an  open  region,  meaning  that  surface  terms 
are  evaluated  as  s  is  approached  from  the  inside  of  the  body.  In  the 
second  of  (29),  scalar  r  denotes  sources  of  thermal  or  chemical  en¬ 
ergy  per  unit  mass,  and  q  g  TxB  is  the  heat  flux  vector  that  is  con¬ 
tinuous  across  all  interfaces  such  that  ([q],n)  =  0  along  s  [58].  The 
total  internal  energy  is  defined  by 

$=  f  pedv  + ^  f  e*edv.  (30) 

Jv  2  Jy 

The  first  term  on  the  right  accounts  for  the  stored  internal  energy  of 
the  body,  denoted  locally  per  unit  mass  by  e.  The  second  term  on  the 
right  side  of  (30)  represents  the  potential  energy  of  the  electric  field 
that  permeates  the  body  and  underlying  vacuum  (i.e.  the  aether). 
The  combined  electromechanical  rate  of  working  is 

jf  t*vds  +  J  (b  +  b)*vdv  +  ^  jTcd^ds 

+  2  [  ptj)dv  +  f  Qdv.  (31) 

dt  Jy  Jy 

The  first  term  on  the  right  of  (31)  accounts  for  the  mechanical  trac¬ 
tion,  the  second  term  accounts  for  body  forces,  the  third  term  ac¬ 
counts  for  the  work  done  by  surface  charges,  and  the  fourth  for 
the  work  of  volumetric  charges.  The  final  term  is  chosen  such  that 
throughout  E3  (moving  dielectric  and  vacuum)  the  balance  of  energy 
is  satisfied  identically: 

a=  [<j),a  -  <j),bVaVb  +  (j):aVbVb]da  +  y  V  bVb ,  (32) 


pe  =  (Vf,Tah  +  ba  -  pva)va  +  <TabVf ,va  -  V; ,qb  +  pr  +  eap  .  (34) 

Then  after  using  (4)  and  (23),  the  local  balance  of  energy  remains 

pe=  (<r,Lr)  -  (v,q)  +  pr+  <e,p>.  (35) 

Notice  that  (35)  differs  from  the  energy  balance  for  non-polar 
continua  in  two  ways.  Firstly,  the  Cauchy  stress  is  not  necessarily 
symmetric  in  (35),  so  that  the  spin  (skew  part)  of  L  may  contribute 

to  the  stress  power  (g,Lt)  =  oabLab  =  aabvbva-  Since  ta  =  oabnb ,  the 
index  of  stress  associated  with  the  traction  vector  is  conjugate  to 
that  associated  with  the  velocity  va,  in  agreement  with  other  non¬ 
linear  theories  for  dielectric  media  [58]  and  also  in  agreement  with 
models  of  generalized  continua  with  couple  stresses  [45].  Secondly, 
the  final  term  on  the  right  side  of  (35)  is  absent  in  non-polar  solids. 

Two  approximations  are  often  made  to  simplify  the  governing 
equations  of  deformable  dielectrics.  The  first  is  the  assumption  of  ge¬ 
ometric  linearity,  i.e.  small  deformations.  In  that  case  introduction  of 
(13)— (1 9)  is  unnecessary,  since  the  distinction  between  undeformed 
and  deformed  configurations  is  not  made  explicitly.  Balances  (23) 
and  (27)  are  unchanged,  but  (35)  becomes  pe  =  <7abVbua  -  Vbqb  + 
pr  +  eap  ,  where  u  is  the  displacement.  The  second  is  the  assump¬ 
tion  that  terms  on  the  order  of  the  product  of  the  electric  field  and 
polarization,  or  on  the  order  of  the  square  of  the  electric  field,  may 
be  neglected  in  the  governing  equations  [54].  In  that  case  body  force 
(20)  and  Maxwell  stress  (21)  vanish,  and  momentum  balances  re¬ 
duce  to  those  of  classical  continua.  With  these  reductions,  the  term 
in  parentheses  in  (34)  still  vanishes,  and  the  final  form  of  the  energy 
balance  (35),  remains  unchanged.  However,  the  stress  tensor  is  now 
symmetric,  and  hence  the  skew  part  of  the  velocity  gradient  does 
not  contribute  to  the  rate  of  change  of  internal  energy. 

Different  contributions  to  the  energy  balance  from  elec¬ 
tromechanical  non-linearity  are  suggested  in  different  theories 
[27,50,58,59];  the  local  balance  (35)  matches  that  of  Tiersten  [58] 
except  for  Tiersten’s  use  of  polarization  per  unit  mass  rather  than 
the  polarization  per  unit  volume.  Polarization  gradients  may  be  im¬ 
portant  for  describing  some  physical  phenomena  [10,49,52];  in  such 
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cases,  augmentation  of  the  energy  balance  to  account  for  effects  of 
polarization  gradients  may  be  necessary. 

2.5.  Entropy  production 


with  reference  heat  flux  Q  =  JF_1q,  second  Piola-Kirchhoff 
stress— possibly  non-symmetric  with  components  IAB  =  JF~1Aoab 

Q  g 

F“1B— and  reference  temperature  gradient  vA0  =  F^va6  all  remain¬ 
ing  unchanged  under  rigid  rotations  [17]. 


The  global  form  of  the  Clausius-Duhem  inequality  is  written 
[45,50,58] 

|36) 

where  p  is  the  entropy  per  unit  mass  and  0(X,  t)  is  the  tempera¬ 
ture.  Application  of  the  divergence  theorem  and  differentiation  of 
the  Helmholtz  free  energy  ^  =  e  -  Op  provides  a  local  form  of  (36) 


3.2.  Thermodynamics 

The  stress  power  entering  (35)  and  (38)  can  be  written  as 

oabLa  b  =  r 1 U FfAacbgac)FaA  =]~ 1  P?F%  (42) 

where  P  is  the  first  Piola-Kirchhoff  stress.  Expanding  the  rate  of  free 
energy  using  (39)  gives 


p{e  -  ij/  -  9t]  -  r)  +  (v,q)  -  l(v0,q)^O.  (37) 

The  entropy  production  inequality  following  from  insertion  of  (35) 
into  (37)  is 

(«.D>  -  (ff.W)  +  <e,p) -p(ils  +  6rj)-  l(v0,q)^O,  (38) 

where  the  covariant  velocity  gradient  L  is  decomposed  into  a  sym¬ 
metric  part  D  and  skew  part  W. 


=  /M,e\  +  /#p\+M 


dE 


\dP 


.P  + 


do 


(43) 


where 


(44, 

dEAB  dEAB  dPA  dPA  dPA 

Substitution  of  (42)-(44)  into  dissipation  inequality  (38)  then  leads 
to 


3.  Elastic  dielectric  solids 

Elastic  dielectrics  considered  in  the  present  section  contain  no 
defects.  Such  materials  obey  the  Cauchy-Born  hypothesis  [6,26],  that 
is,  both  the  material  and  the  primitive  Bravais  lattice  vectors  of  the 
crystal  structure  deform  via  the  deformation  gradient  F  of  (1 ).  Under 
a  homogeneous  deformation  in  the  sense  of  Born  and  Huang  [6], 
a  polarized  dielectric  may  also  exhibit  a  relative  translation  among 
different  atomic  species  (e.g.,  positive  and  negative  ions)  comprising 
the  basis  of  the  crystal  structure.  The  current  Section  provides  a  point 
of  reference  for  comparison  with  the  theory  for  defective  crystals 
developed  later  in  Section  4. 

3.1.  Constitutive  assumptions 


r’p-pgF^-pPF-1®^, 

d  E  dP 


-|,F  +  e-pF-X  p  )-p  U£+t/  0 


8P 


86 


-r’r1  (v6,t 


q)>  o. 


(45) 


Consider  first  isothermal  conditions  for  which  temperature  rates 
and  temperature  gradients  vanish.  Under  such  conditions,  presum¬ 
ing  that  the  dissipation  must  remain  non-negative  when  either  the 
rate  of  spatial  polarization  or  the  rate  of  deformation  gradient  as¬ 
sumes  an  arbitrary  value,  the  following  constitutive  relationships 
follow  from  (45): 


t  r 

E=cPik’ 

8P 


Ea  =  CabP 


8ip_ 

8Pb’ 


(46) 


Constitutive  functions  are  first  assumed  to  exhibit  the  following 
dependencies,  prior  to  consideration  of  objectivity  requirements: 

^  =  ^(F,p,  0),  p  =  p(F,p,  0),  <j  =  <f(F,  p,  0),  e  =  e(F,p,  0), 

q  =  q(F,  p,  0,  v  0).  (39) 

Use  of  polarization  as  an  independent  state  variable  and  electric  field 
as  a  dependent  variable  follows  general  schemes  of  Devonshire  [20] 
and  Toupin  [59].  The  choice  of  electric  field  as  independent  variable 
and  polarization  as  dependent  variable  is  also  possible  [22,50,60]. 
Via  (6),  the  electric  displacement  d  could  substitute  for  either  of  the 
electric  field  or  polarization  in  the  thermodynamic  potentials. 

Consider  rigid  body  motions  of  the  form  x  Qx  +  c,  where 
4  =  0rT  is  a  constant  rotation  tensor  and  c  is  a  constant  translation 
vector.  Under  such  motions,  spatial  polarization  and  electric  field 
vectors  transform  as  p  ->  QP  and  e  ->  Qe,  and  the  remaining  non¬ 
scalar  variables  in  (39)  transform  according  to  F  ->  QF.  Q^Q7. 
q  ->  Qq.  and  0,a  ->  Qab0’b.  On  the  other  hand,  the  referential  po¬ 
larization  and  electric  field  vectors  are  invariant  under  rigid  body 
motions 

P  =  FTp^FTQrQp  =  P,  E  =  FTe^FTQ7Qe  =  E,  (40) 

and  thus  are  valid  candidates  for  use  in  frame-indifferent  constitutive 
relations.  The  following  objective  forms  of  (39)  are  suggested: 

<I/  =  <I/(E,P,6),  ri  =  ri( E,P,0)  £  =  £(E,P,0),  E  =  E(E,P,0), 

Q.  =  Q.(E,  P,  6,  v  9),  (41) 


P=„„  (p|+PF-.  ®  .  P“.p„  .  (47) 

Substituting  (46)  into  (47),  second  Piola-Kirchhoff  and  Cauchy 
stresses  become,  respectively, 

ZAB  =  F-uP°b  =  Po  —  +  ;C“1/1cPcC“1bdEd. 
cEab 

ffab  =J-'FaAZABFbB  =  2p2^-  +paeb,  (48) 

dgab 

where  c~1AB  =  F~ugabF^B.  The  second  terms  on  each  of  the  right 
sides  of  (48)  account  for  possibly  non-symmetric  parts  of  the  stress 
tensors.  These  terms  arise  from  second-order  electromechanical  in¬ 
teractions.  Consistency  of  skew-symmetric  parts  of  (27)  and  (48)  is 
revealed  by  cr^  =  =  T^ba\  and  the  symmetric  total  stress  is 

x  =  2p^  +  p®e  +  e(g)p  +  eoe(g)e  -  ^(e-e)g-1.  (49) 

dg  2 

The  first  term  on  the  right  of  (49)  is  recognizable  as  the 
Doyle-Ericksen  formula  for  hyperelastic  solids  [23,47].  The  final  two 
terms  in  (49)  contribute  to  electrostriction,  even  in  non-polar  solids. 
Presuming  that  (45)-(47)  must  hold  in  the  absence  of  temperature 
gradients 


(50) 
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The  entropy  inequality  thus  reduces  to  the  heat  conduction  inequal¬ 
ity.  When  a  referential  version  of  Fourier’s  Law  applies,  then 


Q  =  -  Kv0,  QA  =  -KAB0<Bt 


(51) 


where  the  conductivity  K  is  symmetric  and  positive  semi-definite 
such  that  -Qa6a  =  Kab0j\6$  ^  0. 

The  energy  balance  (35)  is  now  revisited.  The  specific  heat  ca¬ 
pacity  c  is  introduced,  from  (50)  satisfying 


c_  de  _  de  drj  _  ^d2 i// 
dO  drj  dO  d62 


(52) 


Using  (42)-(51 )  and  multiplying  spatial  energy  balance  (35)  by J  gives 


-eft[Po^  )=PO0'7  =  <V,KV0}  +  Por. 


'89 


(53) 


Carrying  out  the  time  derivative  using  (43)  and  substituting  with 
(52)  results  in 


dil/ 


’  de 


dii/ 


~^dt  »  Po^  I  —  -p°°  I  "aa  )  —  +  0(P*E)  +  0(/,  P 


de 


(54) 


where  the  cross-derivatives  account  for  thermoelastic  and  thermo¬ 
electric  coupling,  respectively, 


P  —  ~Po 


d2i// 

dOdE’ 


1  = 


d2i// 

Po89d  P 


(55) 


modification  to  accommodate  more  general  non-linear  behavior.  For 
example,  a  coupled  dependency  of  dielectric  susceptibility  on  tem¬ 
perature  and  even  powers  of  polarization  is  useful  for  describing 
energy  wells  and  phase  transitions  in  ferroelectric  crystals  [20]  and 
variations  in  the  dielectric  constant  with  applied  voltage  [39].  Al¬ 
lowance  of  dielectric  and  piezoelectric  coefficients  to  vary  with  po¬ 
larization  enables  description  of  hysteresis  under  cyclic  electric  fields 
[64].  A  term  quadratic  in  polarization  and  linear  in  strain  can  be 
added  to  (57)  to  account  for  additional  electrostriction  [49]. 

The  piezoelectric  effect,  Aabc  0,  is  present  in  20  of  the  21  non- 
centrosymmetric  crystal  classes  [19,43].  Of  these,  10  crystal  classes 
possess  a  unique  polar  axis  and  may  exhibit  polarization  in  the  ab¬ 
sence  of  applied  electric  fields  and  applied  strains.  In  the  context 
of  (55),  such  materials  possess  a  non-zero  pyroelectric  coefficient  /, 
leading  to  polarization  with  temperature  variation. 

4.  Dielectrics  with  mobile  vacancies 

In  what  follows,  the  theory  of  non-linear  elastic  dielectrics  of  Sec¬ 
tion  3  is  extended  to  account  for  vacancy  defects,  possibly  mobile 
and  possibly  charged.  Such  a  theory  may  be  used  to  describe  dielec¬ 
tric  thin  film  devices  containing  mobile  oxygen  vacancies,  for  exam¬ 
ple  [14-16,53]  or  silicon  carbide  with  mobile  C  or  Si  vacancies  [3], 
though  it  applies  to  more  general  scenarios  as  well  (e.g.  ceramics 
with  point  defects). 

4.1.  Governing  equations 


Equating  (53)  and  (54),  a  rate  equation  for  the  temperature  emerges 

p0c6  =  (v,  K  v  0)  —  0[  (P,  E)  +  (/,  £)  ]  +  p0r,  (56) 

with  the  first  term  on  the  right  capturing  heat  conduction,  the  second 
and  third  terms  capturing  thermomechanical  and  thermoelectrical 
couplings,  respectively,  and  the  final  term  capturing  non-mechanical 
sources  of  heat  energy. 


3.3.  Representative  free  energy 


Free  energy  function  iff  in  the  first  of  (41 )  is  examined  in  more 
detail  for  illustrative  purposes.  Consider  the  following  form  of  the 
free  energy  per  unit  reference  volume: 


Pol//  =  iE  :  C  :  E  +  1(P,  AP)  +  (P,  A  :  E>  -  (0  -  0O)P  :  E 
-(0-0i)<x,p>  +  r, 
where  the  constant  coefficients  are 


C  =  pQ 


d2  iA 


dEdE 


P  =  -Po 


E=0  .  A  =  p0 
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dOdE 
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l  =  -Po 
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d6dP 


A  =  Po 
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dPdE 


E=0 
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E=0 
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(57) 


(58) 


Here,  CABCD  are  elastic  moduli,  Aab  are  inverse  dielectric  susceptibil¬ 
ities,  Aabc  are  piezoelectric  coefficients,  and  Y  =  Y(0)  is  the  thermal 
energy.  The  following  symmetries  emerge  from  (58): 


0/\BCD  _  qCDAB  _  qBACD  _  qABDC  ^AB  _  ^BA  pAB  _  ^ BA 

Aabc  =  Aacb,  (59) 


meaning  that  the  elastic  moduli  contain  up  to  21  independent  co¬ 
efficients,  the  dielectric  susceptibilities  and  thermal  stress  parame¬ 
ters  up  to  six  independent  coefficients,  and  the  piezoelectric  moduli 
up  to  18  independent  coefficients.  The  coefficients  in  (58)  require 


Maxwell’s  equations  of  electrostatics  and  balances  of  linear  and 
angular  momentum  apply  here,  specifically  local  forms  (10),  (11), 
(23)  and  (27).  However,  the  free  charge  density  (5)  is  extended  to 
delineate  non-vanishing  charges  carried  by  vacancies  in  ionic  solids 
from  other  charges,  and  the  balance  of  energy  and  dissipation  in¬ 
equality  are  modified  to  account  for  electrochemical  energy  supplied 
by  fluxes  and  rates  of  mobile  vacancies.  The  global  rate  of  energy 
exchange  from  thermochemical  sources  in  (29)  is  generalized  to  in¬ 
clude  energy  supplied  by  the  vacancy  flux 

1=  j  prdv-  f(q_,n)  ds  -  f  m~(^~,n)ds,  (60) 

Jv  Js  Js 

where  £  e  TXB  is  the  flux  of  vacancies,  with  dimensions  of  velocity 
per  unit  spatial  volume,  and  m  is  the  scalar  chemical  potential,2 
with  dimensions  of  energy.  Equivalently,  £  may  be  viewed  as  the 
number  of  vacancies  traversing  an  oriented  area  element  nds  per 
unit  time  and  m  the  (electro)chemical  energy  carried  per  (charged) 
vacancy.  The  chemical  potential  can  be  associated  with  the  energy 
required  to  move  an  atom  from  an  interior  lattice  site  to  a  lattice 
site  on  the  surface,  leaving  behind  a  vacancy  in  the  interior  [34].  The 
vacancy  flux  is  continuous  across  interfaces  such  that  ([£],n)  =  0 
along  s.  The  sign  convention  follows  that  used  for  the  heat  flux  q, 
meaning  that  when  m  is  positive,  vacancies  and  energy  flow  out  of 
the  body  when  n)  >0  on  s.  When  the  contacting  medium 

along  s  is  impermeable  to  vacancies,  (£-,  n)  =  (£+,  n)  =  0.  On  the  other 
hand,  when  the  contacting  medium  along  s  is  vacuum  as  opposed 
to  another  solid  body,  vacancies  are  annihilated  as  they  flow  from 
the  dielectric  body  into  the  vacuum  or  are  created  as  they  flow 
from  the  vacuum  into  the  dielectric  body.  It  is  important  to  note 
that  m  is  not  a  constant,  and  instead  may  generally  depend  upon 
field  variables  such  as  stress,  temperature,  and  charge  density.  For 
simplicity,  only  a  single  species  of  vacant  atoms  or  ions  is  considered, 
meaning  that  each  vacancy  carries  the  same  charge  and  occupies  the 


2  More  generally,  a  tensor  chemical  potential  with  components  mab  could  be 
used  [7,31]  so  that  the  energy  flux  becomes  mab(fana.  Here  (60)  is  limited  to  isotropy, 
implying  that  mab  =  mdab. 


682 


J.D.  Clayton  /  International  Journal  of  Non-Linear  Mechanics  44  (2009)  675  -  688 


same  volume.  Mass  transport  associated  with  movement  of  atoms 
between  interstitial  sites  is  not  considered. 

Let  £  denote  the  number  of  vacancies  per  unit  spatial  volume. 
The  vacancy  content  of  the  body  increases  (decreases)  only  when 
vacancies  enter  (exit)  at  interfaces  with  other  bodies  or  free  surfaces. 
Internal  generation  and  annihilation  of  vacancies  are  not  considered. 
Hence,  the  global  time  rate  of  vacancies  per  unit  volume  in  the  body 
in  the  current  configuration  is  dictated  by  the  conservation  law 

jt!fdv=-Sf'n)ds-  (6l) 

Relation  (61)  is  mapped  to  the  reference  configuration  as 

J  Z0dV  =  -  f  «£,N)dS,  (62) 

where  £0  =  ft;  and  =  JF_1£.  Applying  Gauss’s  theorem  to  (62)  and 

localizing  the  result  gives 

to  =  -(V,£0),  (63) 

which  enforces  conservation  of  vacancies  in  the  bulk  solid  [14,15,32]. 

It  is  assumed  that  electric  charges  may  be  carried  by  vacancies 
such  that  their  contribution  to  the  charge  density  in  (5)  is  [14,15,36] 

PV  =  ez£,  (64) 
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Fig.  2.  Kinematic  description  of  deforming  crystal  volume  element  of  fixed  mass; 
dark  circles  represent  atoms  and  open  circles  represent  vacancies. 


Applying  Piola’s  identity  and  (63)  to  the  first  integrand  on  the 
right  of  (67)  results  in 


where  z  is  an  integer  valence  number.  For  vacant  anion  sites,  z  >  0, 
while  for  vacant  cation  sites,  z  <  0.  The  charge  density  (5)  is  modified 
as 

p  =  pc  +  pv  =  ^2  hwezw  +  ez£,  (65) 

i 

where  pc  is  the  density  of  electronic  charges  introduced  in  (5)  and 
pv  accounts  for  the  contribution  of  charged  defects  whose  density 
may  evolve  with  time  via  (61).  Global  charge  conservation  suggests 
that  electronic  charges  contributing  to  pc  in  (65)  could  vary  with 
time  in  order  to  compensate  for  a  non-zero  rate  of  pv.  However,  be¬ 
cause  electronic  charges  travel  much  faster  than  diffusing  vacancies, 
and  because  time  scales  of  present  interest  correspond  to  those  asso¬ 
ciated  with  vacancy  diffusion,  it  is  assumed  that  contributions  from 
moving  electronic  charges  (i.e.,  electric  current)  can  be  neglected  in 
the  governing  equations,  consistent  with  the  quasi-electrostatic  ap¬ 
proximation.  In  other  words,  free  charges  are  assumed  to  rapidly 
self-equilibrate  [62].  With  these  assumptions,  (63)  is  analogous  to 
the  proportionality  relationship  between  the  rate  of  free  electronic 
charges  and  the  negative  divergence  of  an  unsteady  electric  current 
in  electrodynamics  [49].  Note  that  pv  contributes  to  (11)  and  (20). 
Because  the  time  rate  of  pv  may  be  non-zero,  the  local  rate  of  elec¬ 
tromechanical  energy  in  (32)  becomes 

0=  [fa  ~  fbvavb  +  favbvb]da  +  y ^’<‘vbVb 

-4>[f)v +  pvVbvb],  (66) 


-m(v,Q  =  -J  'm(v,t;0)-m(v  '(J  ’F K0)  =J  ’m£o 

=  m(t  +  £trL).  (68) 

Upon  replacing  the  second  of  (29)  with  (60)  and  using  (66)  in 
place  of  (32),  global  energy  balance  (28)  is  consulted.  Then,  after  ap¬ 
pealing  to  (64),  (67)  and  (68),  the  local  balance  of  energy  is  expressed 
as 

pe  =  (<j,  Lr)  -  (v,  q)  +  pr  +  (e,  p)  +  m( £  +  £  trL)  -  <v  m,  Q 

-  ez(j)(t  +  c^trL),  (69) 

differing  from  (35)  via  the  addition  of  the  final  three  terms.  These 
terms  are  associated  with  the  chemical  potential  energy  supplied  by 
vacancy  rates  and  fluxes  and  the  electromechanical  energy  attributed 
to  the  time  rate  of  change  of  charged  vacancies.  Substituting  (69)  into 
(37)  provides  the  dissipation  inequality  for  dielectrics  with  charged 
vacancies,  referred  to  configuration  B 

(a,  Lt)  +  (e,  p)  +  m{t  +  £  trL)  -  (v  m,  0  -  ez0(^  +  £  trL)  -  p( \j/  +  6rj) 
-l(V0,q)^O.  (70) 

4.2.  Kinematics 

The  deformation  gradient  of  (1)  is  decomposed  multiplicatively 
[4,44]  to  account  for  recoverable  thermoelastic  deformation  FE  and 
volumetric  deformation  Fv  attributed  to  vacancies 

F  =  FeFv,  FaA  =  FE“Fy«,  (71) 


where  the  final  term  on  the  right  side  accounts  for  the  increase 
in  total  internal  energy  resulting  from  an  evolving  charge  density 
associated  with  mobile  vacancies. 

From  the  divergence  theorem,  the  chemical  work  associated  with 
the  final  term  on  the  right  of  (60)  is 


n)  ds  = 


(67) 


implying  the  existence  of  an  unstressed  intermediate  configuration 
[25]  labeled  in  Fig.  2  as  B,  and  a  corresponding  series  of  global  tan¬ 
gent  mappings  Fv  :  TB0  -*  TB  and  FE  :  TB  -*  TB  [4,13].  The  elastic 
deformation  FE  consists  of  stretch  of  the  lattice  due  to  mechanical 
loading,  elastic  and  rigid  body  rotations,  and  deformation  result¬ 
ing  from  thermal  expansion  or  contraction  induced  by  temperature 
changes.  The  latter  is  considered  reversible:  in  the  absence  of  other 
deformation  modes,  restoration  of  the  local  crystal  element  to  its 
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reference  temperature  relieves  the  thermal  strain.  Deformation  from 
vacancies  Fv  is  mechanically  irreversible,  since  it  generally  remains 
when  external  forces  are  removed  from  the  crystal,  and  is  assumed 
to  manifest  only  through  a  change  in  volume  of  the  solid.  Possible 
shape  changes  of  the  solid  associated  with  non-spherical  pores  or 
anisotropic  clusters  of  vacancies  within  a  given  volume  element  of 
crystal  are  not  considered.  Because  the  deformation  resulting  from 
vacancies  is  isotropic,  it  can  be  expressed  as 

FV  =  (1  —  0)-1/36.  =  (1  -  0)“1/34  (72) 

The  scalar 

<t>  =  {dV-  dV)/dV  =  1  -Jv~'  (73) 

is  the  volume  fraction  of  vacancies  per  unit  volume  in  configuration 
B  [2,13,42],  where  Jv  is  the  Jacobian  determinant  of  Fv.  Neither  FE 
nor  Fv  need  be  integrable  to  a  displacement  function;  specifically, 
incompatibility  of  the  latter  is  reflected  by  the  quantity  S^A(j)  B],  non¬ 
zero  when  the  field  <£(X,  t)  is  heterogeneous.  An  external  Cartesian 
coordinate  frame  with  basis  vectors  not  tangent  to  any  material  lines 
is  implied  for  anholonomic  configuration  B  [12],  with  corresponding 
Cartesian  metric  g =  S^.  Also,  5  is  the  shifter  between  B  and  B0 , 
equivalent  in  (73)  to  the  identity  map  1  with  coincident  coordinate 
axes  implied  in  B  and  B0. 

As  illustrated  in  Fig.  2,  the  reference  configuration  is  a  perfect 
crystal,  free  from  defects,  though  the  theory  does  not  preclude  a  dis¬ 
tribution  of  vacancies  at  t= 0.  Notice  that  by  (72),  vacancies  increase 
the  specific  volume  and  decrease  the  mass  density  of  an  element  of 
fixed  mass  when  </>  >  0.  As  shown  in  Fig.  2,  the  total  number  of  atoms 
in  the  element  under  consideration  remains  fixed  in  all  configura¬ 
tions,  i.e.  mass  is  conserved  in  agreement  with  (2).  Introduction  of  a 
vacancy  by  removal  of  an  atom  from  a  perfect  lattice  site  can  cause 
contraction  of  the  remaining  atoms  towards  the  defect  such  that  the 
energy  associated  with  interatomic  separation  distances  is  reduced, 
shown  qualitatively  by  the  local  attraction  of  a  few  of  the  atoms  to 
their  neighboring  vacant  sites  in  Fig.  2.  For  example,  according  to 
continuum  elasticity  theory,  the  radial  displacement  field  induced  by 
a  spherical  point  defect  in  an  infinite  isotropic  body  decays  inversely 
to  the  square  of  the  radial  distance  from  the  defect  [28].  When  the 
reduction  in  volume  induced  by  such  local  displacement  fields  (i.e. 
the  relaxation  volume)  is  small  compared  to  the  volume  associated 
with  each  vacant  lattice  site  (i.e.  the  atomic  volume),  then  Jv  >  1. 
Otherwise,  when  the  relaxation  volume  is  larger  than  the  atomic 
volume  due  to  strong  interatomic  forces,  Jv  <  1  and  (j)  <  0,  and  the 
overall  volume  per  unit  mass  relative  to  a  perfect  crystal  is  reduced 
by  vacancies.  For  ionic  crystals,  additional  expansion  or  contraction 
resulting  from  electrostatic  interactions  between  ions  surrounding 
charged  vacant  sites  is  conceivable;  such  effects  are  embedded  in  Fv. 
Charge  neutrality  need  not  occur  in  each  volume  element  of  fixed 
mass. 

From  (71)  and  (72),  the  deformation  rate  resulting  from  vacancy 
flux  is 

Lv  =  pVpV-i  =  ^(3  _  (74) 

Let  <J)0  =  ag o  denote  the  fraction  of  vacancies  per  unit  reference 
volume,  where  a  is  a  scalar  constant  denoting  the  volume  change  ef¬ 
fected  by  a  single  vacancy.  Net  expansion  is  represented  by  a  >  0,  and 
net  contraction  by  a  <  0.  The  following  relationships  then  emerge 
between  dimensionless  volume  fractions  and  number  densities  of 
vacancies: 

4>  =JV -1  <Po  =/"  Vo  =  (1  -  4>Mo  =JV,  (75) 

where  ]E  is  the  Jacobian  determinant  of  FE. 


The  following  notation  is  introduced,  to  label  variables  couched 
in  intermediate  configuration  B,  the  configuration  acting  here  as  an 
evolving  reference  state  for  the  thermoelastic  response  [11,25,44]: 


P  =  a/-’=p/.  (76) 

2Ee  =  Ce  — 1,  (*,  =  F*gablf,  (77) 

£=/-1FvEFVT=/FE-1ffFE-T,  (78) 

V6  =  (v  0)FV_1  =  (v  0)FE,  Vm  =  (v  m)Fv_1  =  (v  m)FE,  (79) 

q  =/-1  FVQ  =  /Fe_1  q,  C=V  FVC0  =J£FE-'  (80) 

p  =  FETp,  e  =  FE7e,  d=/FE1d.  (81) 


The  mass  density  per  unit  intermediate  volume  in  B  is  written  as 
p  in  (76).  Introduced  in  (77)  is  the  covariant  elastic  strain  Ee.  In 
(78),  relationships  between  contravariant  elastic  stress  £  eTB  x  TB, 
second  Piola-Kirchhoff  stress  £,  and  Cauchy  stress  a  are  given.  None 
of  these  stress  measures  need  be  symmetric  in  polarized  media. 
The  intermediate  temperature  gradient  VO  and  chemical  potential 

gradient  Vm  follow  in  (79),  where  Va  =  VaFEa  =  VaFVoT1A  defines  the 
anholonomic  covariant  derivative  of  a  scalar  function  [13].  In  (80), 
the  heat  flux  and  vacancy  flux  are  each  mapped  from  the  reference  to 
the  intermediate  configuration  via  appropriate  Piola  transformations. 
Polarization,  electric  field,  and  electric  displacement  are  mapped  to 
B  according  to  (81);  substitution  into  (6)  then  provides 

d=/-1CE-1(s0e  +  p).  (82) 

4.3.  Constitutive  assumptions 

Constitutive  functions  used  here  are  of  the  same  form  as  (39), 
apart  from  five  major  differences  that  reflect  the  presence  of  inelastic 
deformation  and  defects.  The  first  is  that  the  recoverable  thermoe¬ 
lastic  strain  Ee,  as  opposed  to  the  total  strain  E,  is  used  as  an  inde¬ 
pendent  state  variable,  following  the  logic  of  multiplicative  plasticity 
theory  [11,44].  As  such,  the  elastic  stress  £,  work  conjugate  to  Ee, 
is  used  in  place  of  the  second  Piola-Kirchhoff  stress  £.  The  second 
is  that  electric  field  and  polarization  defined  with  respect  to  config¬ 
uration  B  via  (81)  are  used,  as  opposed  to  their  referential  counter¬ 
parts  E  and  P.  The  third  is  that  intermediate  temperature  gradient 
V6  of  (79)  and  heat  flux  q  of  (80)  are  used  as  opposed  to  analogous 
quantities  defined  on  B0.  The  fourth  is  that  the  response  functions 
are  assumed  to  depend  upon  the  density  of  defects  and  the  fifth 
is  that  a  constitutive  function  for  the  vacancy  flux  is  posited  that,  in 
addition  to  possibly  depending  on  the  other  independent  state  vari¬ 
ables,  also  depends  upon  the  gradient  of  the  chemical  potential,  Vm. 
Constitutive  functions  are  written  generically  as 

i^=^(Ee,p,^,0),  tj  =  f/(EE,p,  (*,  6),  E  =  £(Ee,p,£,0),  e  =  e(EE,p,  9), 

Q  =  Q(Ee,  p,  9,  V  9),  £  =  £(EE,p,  9,  Vm).  (83) 

Verification  is  straightforward  that  all  variables  in  (83)  remain 
invariant  under  rigid  body  motions  of  the  form  x  Qx  +  c  with 
Fe  ->►  QFe.  The  dependence  on  number  of  defects  per  unit  spatial 
volume  £  could  alternatively  be  replaced  by  a  dependence  on  the 
volume  fraction  of  defects  <fi  via  (75). 

4.4.  Thermodynamics 

Use  of  (71)-(75)  leads  to  the  expansion  of  the  spatial  velocity 
gradient 


L  =  FeFe-’  +G4/3)[£  +  £tr(FEFE-’)]1,  A  =JEa(1  -JEa£)~\  (84) 
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where  the  scalar  A  appears  prominently  in  what  follows.  The  stress  Gab  _jE- i  pEa^PpEb  _  2 +paeb  -^g°b 

power  entering  (69)  and  (70)  can  be  written  ^  dg„b  d( 


(95) 


cabLab  =  (JE-'h* -Ap{FEa-'x)FExa -Apt  (85) 

where  Paa  =  JEFE~lclaab  is  the  elastic  first  Piola-Kirchhoff  stress,  and 
p  =  -ff"/3  is  the  Cauchy  pressure.  Expanding  the  rate  of  free  energy 
from  (83) 


#  ce\  ,  /#  *\  ,  ,  #■ 


/  +  w7  +dee+dil 


where  from  definitions  (77)  and  (81) 


pE  pEb  #  pEa  #  s  UXV  f,  p Ea  ,  ^  pEap 

t~a  =  r  «  T^gab^.a  •  ~^Poi  =  —Parx  +  —  f  a  Pa- 


0E£  a/?  ^  8EE„ 

ap  a  p 


8x1/  „ 
5p« 


3p« 


Notice  also  that  from  (84) 

^rL  =  «l+^)F£-laF£“+A«. 


(86) 


(87) 


(88) 


Substituting  (84)-(88)  into  dissipation  inequality  (70)  then  yields 


/-1  P^-Ap^FE~u-pFEb  fpfgab-p  ^-Pa-ezkiUA^F^ 


+  m«l+A^)F£a 


FEa+  m(l+A£)-ez</)(l+A£)-Ap-p^ 


dH 


+ 


F? 

8pa 


dip 

Pa~ 

r»+,J 

-J£_1  d  1  Va0g“  - J£-’  Vamf  >  0. 


(89) 


Consideration  of  distinct  thermodynamic  processes  in  conjunction 
with  (89)  leads  to  a  number  of  constitutive  laws  [10,17,30,50].  Pre¬ 
suming  that  (89)  must  be  satisfied  for  arbitrary  rates  of  temperature 
leads  to  the  usual  relationship  between  temperature  and  entropy, 
identical  to  (50).  Presuming  that  dissipation  associated  with  the  rate 
of  polarization  must  be  non-negative  for  arbitrary  rates  of  polariza¬ 
tion  leads  to 


ea  =  F£“p 


8iP_ 

dpx' 


IpP 


djj_ 

dp/)’ 


(90) 


Assuming  that  dissipation  associated  with  the  rate  of  vacancy  con¬ 
centration  must  remain  non-negative  for  arbitrary  concentration 
rates  in  the  absence  of  gradients  in  chemical  potential  or  tempera¬ 
ture,  chemical  potential  m  should  satisfy 


m  =  (1  +A£)  V^  +  ez</)+A(l  +  Ad)  V  (91) 

Requiring  non-negative  dissipation  associated  with  the  elastic  de¬ 
formation  gradient  rate  then  gives 


JE~'hX  =  pFEe-Srgab  +Ap^F£-la  +  p^Pa  +  ez<k(  1  +M)FEa~U 

ft  °P* 


-  m£(  1  +A£)FEa  la, 

which  upon  substitution  of  (91)  becomes 

p.a  -  pEb-nd'l'  ,F .E-la*# 

Pa=gabF,P^r+PaPWx-tFM  p-. 

Using  (90),  the  elastic  stress  and  Cauchy  stress  are  then,  respectively, 


(92) 


(93) 


=  Fe~^  =  pJfL.  +  JECE~^pxC 


8x1/ 


8Ee  R 


rE-XpSp 


-  ,fc£_la^h  — 
F  8C 


The  final  terms  on  the  right  of  each  of  (94)  and  (95)  account  for  pos¬ 
sible  effects  of  vacancy  concentration  on  hydrostatic  pressure.  From 
(95),  consistency  of  skew-symmetric  parts  of  Cauchy  and  Maxwell 
stresses  in  (27)  is  revealed  by  cr^  =  e[bpa^  =  $ba\  and  the  total  stress 
is  thus  symmetric 

x  =  2p-C  +  £p-Cg  1  +p<g>e  +  e(g>p  +  £oe(g)e-  —  (e-e)g  7  (96) 
dg  d£  2 

The  remaining  terms  in  dissipation  inequality  (89),  upon  consul¬ 
tation  of  (90)-(92),  are 

-/-10-1(V0,q>  -  J£_1  (Vm,Q  >0.  (97) 

Though  more  general  formulations  are  possible,  the  contribution 
from  the  flux  of  thermal  energy  is  non-negative  when  Fourier-type 
behavior  occurs  in  the  elastically  unloaded  configuration,  analogous 
to  (51) 

q  =  -KV0,  (Vct0)K^(V/!0)^O,  (98) 

where  thermal  conductivity  matrix  K  e  TB  x  TB  is  symmetric  and 
positive  semi-definite.  Analogously,  the  contribution  to  entropy  pro¬ 
duction  from  the  flux  of  vacancies  can  be  made  non-negative  by 
assuming 

\  =  -DVm,  (Vam)Da^(V|gm)^0,  (99) 

with  diffusivity  matrix  D  e  TB  x  TB  symmetric  and  positive  semi- 
definite.  Applying  (98)  and  (99)  in  (97),  the  dissipation  inequality  is 
always  satisfied: 

(/£0)_1(V0,KV0)  +/“1(Vm,DVm)^0.  (100) 


Finally,  from  the  form  of  the  chemical  potential  in  (91),  the  va¬ 
cancy  flux  in  (99)  satisfies 


r=  -  DaP 


(i  +M)~'p% 


+  ez' 


+v£ 


[A(i+AH)-'p]  J. 


(101) 


The  first  term  on  the  right  of  (101)  represents  effects  of  concentra¬ 
tion  gradients  on  the  flux  of  neutral  or  charged  vacancies,  the  sec¬ 
ond  represents  effects  of  electrostatic  potential  gradients  on  the  flux 
of  charged  vacancies,  and  the  third  term  represents  effects  of  pres¬ 
sure  gradients  on  the  flux  of  neutral  or  charged  vacancies.  Effects  of 
gradients  of  concentration  and  pressure  on  diffusion  are  well  docu¬ 
mented  [33,34,38],  as  are  effects  of  gradients  of  electrostatic  poten¬ 
tial  in  ionic  crystals  [18,37,40,63]. 


4.5.  Representative  free  energy 


Further  insight  is  achieved  by  considering  a  free  energy  per  unit 
intermediate  volume  of  the  form 


pi]/  —  2®^'  ^  •  EE  +  2  (P*  Ap)  H-  (p,  A  :  Ee) 

+  IH2  -  (9  -  0O)P  :  EE  -  (0  -  0,  )(x.p>  -  pc0 In  (Ji)  ,  (102) 


where  Oo  and  0\  are  constants  with  dimensions  of  temperature. 
The  remaining  coefficients  in  (102)  are  isothermal  constants  defined 
according  to 


C  =  p 


d2i]/ 

<3EE  <3EE 


ee=o 

p=o 

£=o 

e=o0 


,  .  aV 
A-',as# 


Ee=0 

p=o 

<^=o 

0=0q 


,  -  02I A 

P0p0EE 


Ee=0  ’ 
p=o 

<^=o 

6=60 


(94) 
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_ 

D  =  p - y 


Ee=0  - 

p=o 

£=o 


c  =  -0 


gV 

302 


Ee=0 

p=0 

^=o 


p=-p 


<30<3EE 


Ee=0  « 

p=o 

0=00 


_  d2i// 
7  =  — P - — 


ee=o  ■ 

p=0 

£=o 

0=0! 


(103) 


is  often  referred  to  as  the  dielectric  constant  [39],  its  entries  here  in 
the  context  of  geometric  non-linearity  are  not  fixed  constants,  but 
instead  depend  upon  the  elastic  deformation.  The  materially  linear 
framework  of  (102)-(108)  is  insufficient  to  address  ferroelectric  be¬ 
havior,  which  requires  a  more  detailed  energy  functional  to  account 
for  transition  temperatures  and  higher-order  influences  of  polariza¬ 
tion  [20]. 


Terms  in  (102)  referred  to  the  intermediate  configuration  are  anal¬ 
ogous  to  those  of  the  theory  for  elastic  dielectrics  in  (57)  referred  to 
the  reference  configuration,  apart  from  the  quadratic  form  of  the  va¬ 
cancy  energy,  D£2/2,  that  is  newly  introduced  in  (102).  Consequences 
of  this  particular  choice  of  vacancy  energy  with  regards  to  diffusion 
are  examined  in  more  detail  in  Section  4.6.  Arguments  regarding 
symmetry  analogous  to  (59)  apply  here  as  well. 

From  (94)  and  (102),  the  elastic  stress  is 

£*  =  +  h-Ax«P  _{e_  _  c£_ ,a^2 

+JECE~uxpxCE-'lPsed,  (104) 

where  the  first  four  terms  on  the  right  side  are  symmetric  contri¬ 
butions.  The  first  term  on  the  right  side  accounts  for  hyperelastic¬ 
ity,  the  second  term  accounts  for  the  inverse  piezoelectric  effect,  the 
third  accounts  for  the  thermoelastic  effect  (i.e.  thermal  expansion  or 
contraction),  the  fourth  accounts  for  vacancies,  and  the  final  term 
reflects  the  possibly  non-symmetric  contribution  of  the  Maxwell 
stress.  When  the  vacancy  concentration  is  small  and  £2  0,  the 

contribution  of  vacancies  to  the  stress  from  the  fourth  term  on  the 
right  of  (104)  is  negligible.  On  the  other  hand,  consider  cases  when 
this  term  is  significant,  when  v  >  0,  and  when  the  elastic  strain  is 
held  fixed.  Under  such  constraints,  the  concentration  of  vacancies 
induces  a  negative  hydrostatic  component  of  £,  or  a  positive  hydro¬ 
static  pressure. 

From  (90)  and  (102),  the  electric  field  referred  to  the  intermediate 
configuration  is 

=f-'(%\Atopx  +  2*%  -  {9 -  8, )*/>].  (105) 

The  first  term  in  braces  on  the  right  side  of  (105)  accounts  for  the  di¬ 
electric  permittivity,  the  second  accounts  for  the  piezoelectric  effect, 
and  the  third  term  accounts  for  the  pyroelectric  effect.  From  (82) 

¥  =  (80 Afa  +/c£-1to)pz  +  -{8-9 1  )xp.  (106) 

Assuming  that  A is  positive  definite  and  inverting  (105) 

px  =  (JEA~)cE~^)ea  -  U~)AfiSe)EEe  +  {8-8,  Ufpf) 

=  2z-“ea- V84  +  (0-0,)ir  (107) 

where  A  depends  on  the  elastic  strain  and  where  A  and  7  are  con¬ 
stants.  Substitution  of  (107)  into  (106)  then  yields  the  relationship 
between  electric  displacement  and  electric  field 

¥  =  [eaJECE-'P“  +  JE2CE-^A~JcE-u«\ea  -  [(s0  -  1  )Apde 

+]ece-^a;Ia^\ees  +  k  80  -  +JECE-^A~Jf](e-  9,) 

=  Eo4*ea-Ali5eEE£  +  xIS(6-e1),  (108) 

where  sf'  =  jece-]7(’  +  epJE2CE^]7/A  f  C£~ 1  ‘A1  is  the  relative  per- 
mittivity  and  is  clearly  symmetric,  and  where  the  other  notation  is 
evident  from  (108).  The  first  term  on  the  right  of  (108),  £o£rC,  rep¬ 
resents  the  purely  dielectric  effect,  the  second  A  :  Ee  accounts  for 
piezoelectric  coupling,  and  the  third  %(0  -  6 1)  accounts  for  pyro¬ 
electric  coupling.  Regarding  the  latter,  temperature-induced  polar¬ 
ization  will  occur  when  7  is  non-zero  and  6±6\  in  (107).  While  8r 


4.6.  Diffusion  of  vacancies 

Next  consider  diffusion  of  vacancies  in  the  context  of  chemical 
potential  (91),  diffusion  law  (101),  and  free  energy  function  (102). 
The  chemical  potential  m  in  this  case  is 

m  =/-1(l  +  A^)-1^  +  ez0  +/ap,  (109) 

leading  to  the  flux  equation 

f  =  -D«h\JE~\  1  +AO-'m,p  +  ezfli+JEocpflj 

+  t\JE-\  1  +Ai)-1b],p  +p\JEa]p}.  (110) 

The  first  term  on  the  right  of  (110)  causes  vacancies  to  diffuse  from 
regions  of  high  concentration  to  regions  of  low  concentration.  The 
second  term  causes  positively  charged  vacancies  to  diffuse  from  re¬ 
gions  of  high  electrostatic  potential  to  regions  of  low  electrostatic 
potential,  leading  to  a  reduction  in  system  energy  since  the  vacancy 
charge  density  ez£  and  electrostatic  potential  </>  are  work  conju¬ 
gates  in  (66).  The  third  term  causes  vacancies  to  diffuse  from  regions 
of  high  hydrostatic  pressure  to  regions  of  low  pressure  a  >  0.  The 
fourth  and  fifth  terms  reflect  non-linear  effects  resulting  from  spatial 
gradients  in  property  v  (e.g.  heterogeneous  crystals),  elastic  volume 
changes,  and  vacancy  concentration.  Upon  neglecting  non-linear  ef¬ 
fects  such  that  1  +  A£  &  1  in  a  homogeneous  body  in  the  absence 
of  elastic  deformation,  pressure  gradients,  and  electric  fields,  (110) 
reduces  to  a  version  of  Fields  Law 

f«-6  (111) 

i.e.  a  familiar  linear  relationship  between  the  flux  of  vacancies  and 
the  spatial  gradient  in  vacancy  concentration.  The  dependence  of  free 
energy  upon  vacancy  density  used  in  (102)  is  perhaps  the  simplest 
formulation  that  maintains  elements  of  physical  realism.  The  partic¬ 
ular  form  used  for  the  vacancy  energy,  or  the  parameter  D,  can  be  es¬ 
timated  using  arguments  from  chemical  mixture  theory  [29]  or  free 
energy  versus  defect  content  results  from  lattice  statics/dynamics 
calculations  [11],  at  defect  densities  large  enough  such  that  interac¬ 
tion  effects  between  neighboring  defects  become  significant.  More 
general  functions  than  (102)  are  required  to  address  all  relevant  cou¬ 
plings  between  defect  content  and  bulk  electromechanical  and  ther¬ 
mal  responses.  For  example,  the  vacancy  content  in  elastic  dielectrics 
may  alter  the  elastic  moduli  [14,21]  and  affect  the  piezoelectric  co¬ 
efficients.  Furthermore,  the  vacancy  energy  and  hence  diffusion  may 
be  strongly  affected  by  temperature,  and  could  perhaps  more  realis¬ 
tically  be  a  non-quadratic  function  of  temperature  and  vacancy  con¬ 
centration  in  dielectric  solids  [14,15,63].  In  the  present  treatment, 
no  restrictions  are  placed  on  the  diffusion  coefficients  D0^,  except 
that  their  matrix  should  be  symmetric  and  positive  semi-definite  to 
ensure  a  positive  rate  of  entropy  production.  It  is  well  known  that 
entries  of  D ^  need  not  be  constants,  but  instead  may  depend  upon 
pressure  and  temperature  [38].  Deviatoric  stresses  can  affect  the  dif¬ 
fusion  of  point  defects,  as  observed  for  creep  mechanisms  in  poly¬ 
crystals  [33]  and  vacancy  or  interstitial  migration  in  semiconductors 
[1].  Effects  of  stress  directionality  on  diffusivity  could  be  incorpo¬ 
rated  via  specification  of  D dependent  on  the  elastic  strain. 
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Fig.  3.  Dielectric  slab  with  vacancies  subjected  to  biaxial  lattice  strain  and  applied 
electrostatic  potential. 


4.7.  Example:  diffusion  under  biaxial  lattice  strain 

Consider  the  slab  of  material  shown  in  Fig.  3,  subjected  to  the 
uniform  biaxial  elastic  strain  field  sE  in  the  x1  -  x2  plane,  so  that  in 
Cartesian  coordinates,  F^1  =  F|2  =  1  +  sE  and  FEf  =  Saa  otherwise.  Po¬ 
tentials  are  applied  to  faces  of  the  slab  normal  to  x3,  with  values  of 
</>  differing  on  either  face.  Thus  the  voltages  impose  an  electric  field 
-V0,  where  V  is  the  spatial  gradient  operator  in  the  x3 -direction. 
For  example,  when  one  face  of  the  slab  is  grounded  with  cj)  =  0, 
the  magnitude  of  the  applied  electric  field  is  then  the  value  of  ^ 
at  the  opposite  face  divided  by  the  thickness  of  the  slab  in  the  x3- 
direction.  For  simplicity,  the  interfaces  here  are  assumed  to  prevent 
global  straining  in  the  x3 -direction.  The  system  corresponds  physi¬ 
cally  to  a  dielectric  film  device  [14-16,53].  In  the  physical  system, 
the  elastic  strain  sE  may  arise  from  a  mismatch  in  lattice  parameters 
across  the  interfaces,  for  example  between  the  dielectric  slab  and 
neighboring  substrates  or  electrodes,  in  conjunction  with  processing 
steps  that  may  introduce  residual  thermal  stresses  [14,53].  Consider 
uniaxial  diffusion  with  vacancy  flux  (  parallel  to  the  x3 -direction,  in 
a  dielectric  solid  with  isotropic  diffusivity  d  such  that  D0^  =  dS 
Under  these  conditions,  flux  equation  (110)  mapped  to  the  spatial 
coordinate  frame  becomes 

-C/d  =  [(1  +  sEr2  -  2oc{JWt  +  [(1  +  sEr2]ezV(j),  (112) 

where  o=JE-^o=pd2\l//d£2.  Terms  in  braces  in  (112)  reduce  to  unity 
in  the  context  of  traditional,  geometrically  linear  models  in  which 
finite  deformations  are  not  considered  [18,37],  wherein  the  chemical 
potential  often  exhibits  the  simple  form  m  =  v£  +  ezcf  +  a p.  In  the 
scenario  corresponding  to  Fig.  3,  a  spatially  constant  pressure  p  is 
assumed  to  complement  the  assumption  of  a  constant  elastic  strain; 
hence  possible  effects  of  pressure  gradients  on  diffusion  are  absent 
in  (112). 

Shown  graphically  in  Fig.  4(a)  is  the  normalized  flux  ~C/(dvV£) 
corresponding  to  the  first  term  on  the  right  of  (112).  Compres¬ 
sive  strains  sE  <  0  and  reductions  in  volume  from  vacancies  <  0 
tend  to  increase  the  normalized  flux,  while  tensile  strains  sE  >  0 
and  positive  volume  changes  a£>0  from  vacancies  tend  to  de¬ 
crease  the  normalized  flux.  Interestingly,  at  very  large  concentra¬ 
tions  a^0.2  and  tensile  strains  SE>  0.6,  the  instantaneous  direc¬ 
tion  of  flux  would  change  since  then  -f/(dbV^)  <  0,  and  vacancies 
would  coalesce  rather  than  migrate  to  areas  of  lower  concentration. 
Shown  in  Fig.  4(b)  is  the  normalized  flux  arising  only  from  gradients 
in  electrostatic  potential,  corresponding  to  the  second  term  on  the 
right  of  (112).  Again,  compressive  strains  SE<  0  tend  to  increase  the 
flux  while  tensile  strains  sE  >  0  tend  to  decrease  the  flux.  In  sum¬ 
mary,  large  elastic  deformations  \sE\^0A  and  large  vacancy  concen¬ 
trations  |a^|^0.1  are  predicted  according  to  (112)  and  Fig.  4  to  influ¬ 
ence  vacancy  migration.  However,  stress  and  vacancy  content  may 
also  affect  the  diffusivity  d  [38];  such  effects  would  depend  on  the 
particular  dielectric  material  comprising  the  slab  and  are  not  consid¬ 
ered  in  Fig.  4,  wherein  d  is  assumed  constant.  Furthermore,  disloca¬ 
tion  nucleation  or  fracture  could  ensue  in  dielectric  solids  to  relieve 
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Fig.  4.  Normalized  vacancy  flux  driven  by  (a)  concentration  gradient  and  (b)  elec¬ 
trostatic  potential  gradient. 


potentially  large  elastic  strains;  such  effects  are  also  not  considered 
here. 

5.  Conclusions 

A  model  framework  has  been  formulated  for  the  electromechan¬ 
ical  behavior  of  dielectric  crystalline  solids  subjected  to  large  defor¬ 
mations.  These  deformations  consist  of  thermoelastically  recoverable 
stretch  and  rigid  body  rotation,  as  well  as  volumetric  changes  result¬ 
ing  from  voids  or  vacancies.  Charge  transport  resulting  from  diffusion 
is  considered;  the  motion  of  electronic  charges  occurring  at  faster 
time  scale  (i.e.,  electric  currents)  is  not.  A  model  for  elastic  dielectrics 
without  defects,  in  the  quasi-electrostatic  approximation,  has  been 
presented  first,  to  lend  context  to  the  more  general  theory  allow¬ 
ing  for  vacancies  developed  later.  The  latter  theory  may  be  the  first 
to  consider  finite  volume  changes  resulting  from  charged  vacancies 
in  ionic  crystals  in  the  context  of  multiplicative  inelasticity.  Govern¬ 
ing  equations,  constitutive  relations,  and  kinetic  equations  are  cast 
in  the  relaxed  intermediate  configuration  implied  by  the  multiplica¬ 
tive  decomposition.  Non-linear  electromechanical  effects  associated 
with  Maxwell’s  stress  are  retained.  Effects  of  spatial  gradients  of  va¬ 
cancy  concentration,  electrostatic  potential,  and  hydrostatic  pressure 
arise  naturally  in  the  thermodynamically  admissible  diffusion  equa¬ 
tion  for  the  vacancy  flux,  without  assumption  of  such  kinetic  law 
dependencies  a  priori.  Upon  assumption  of  a  quadratic  free  energy 
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dependence  on  vacancy  density,  Fick-type  diffusion  is  recovered  in 
the  linearized  limit,  in  the  absence  of  electrostatic  and  pressure  ef¬ 
fects.  However,  in  the  non-linear  regime,  large  elastic  deformations 
and  large  vacancy  concentrations  may  influence  the  vacancy  flux. 
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valence  of  carrier  i 
no.  of  carriers  per  spatial  volume 
charge  magnitude  of  an  electron 
permittivity  of  vacuum 
free  surface  charge  density 
internal  surface  charge  density 
reference  free  charge  density 
spatial  free  charge  density 
charge  density  of  vacancies 
charge  density  of  electrons/holes 
reference  polarization 
spatial  polarization 
intermediate  polarization 
reference  electric  displacement 
spatial  electric  displacement 
intermediate  electric  displacement 
mechanical  body  force 
electromechanical  body  force 
mechanical  traction 
total  traction 
Cauchy  stress 
Maxwell  stress 
total  stress 

first  Piola-Kirchhoff  (P-I<)  stress 
elastic  first  P-I<  stress 
second  P-I<  stress 

elastic  second  P-I<  stress 
Cauchy  pressure 
kinetic  energy 

total  internal  (system)  energy 
thermochemical  energy 
total  external  power 
electrostatic  power  per  unit  volume 
absolute  temperature 
free  energy  per  unit  mass 
internal  energy  per  unit  mass 
entropy  per  unit  mass 
heat  source  per  unit  mass 
specific  heat  per  unit  mass 
thermal  energy 

reference  vacancy  volume  fraction 
intermediate  vacancy  fraction 
no.  vacancies  per  reference  volume 
no.  vacancies  per  spatial  volume 
net  volume  per  vacancy 
normalized  volume  per  vacancy 
chemical  potential 
reference  vacancy  flux 
spatial  vacancy  flux 
intermediate  vacancy  flux 
intermediate  vacancy  diffusivity 
isotropic  diffusivity 
reference  heat  flux 
spatial  heat  flux 
intermediate  heat  flux 
reference  thermal  conductivity 
intermediate  thermal  conductivity 
spatial  thermal  stress  coefficients 

intermediate  thermal  stress  coefficients 
spatial  pyroelectric  coefficients 
intermediate  pyroelectric  coefficients 
reference  elastic  coefficients 
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C,  C  intermediate  elastic  coefficients 

A,  Aab  reference  dielectric  susceptibility 

A,  intermediate  susceptibility 

A,  Aabc  reference  piezoelectric  coefficients 

A,  A^1  intermediate  piezoelectric  coefficients 

£R,  ejjf  relative  dielectric  permittivity 

v  spatial  vacancy  energy  coefficient 

v  intermediate  vacancy  energy  coefficient 

c  intermediatespecific  heat  coefficient 

0o,  0\  reference  temperatures 
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